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Abstract 

Motivated by the way Japanese tatami mats are placed on the floor, 
we consider domino tilings with a constraint and estimate the number 
of such tilings of plane regions. We map the system onto a monomer- 
dimer model with a novel local interaction on the dual lattice. We make 
use of a variant of the Hamiltonian replica exchange Monte Carlo method 
where data for ferromagnetic and anti-ferromagnetic models are combined 
to make a single family of histograms. The properties of the density of 
states is studied beyond exact enumeration and combinatorial methods. 
The logarithm of the number of the tilings is linear in the boundary length 
of the region for all the regions studied. 


1 Introduction 

A domino tiling is a non-overwrapping covering of a planar region with rect¬ 
angular tiles of edge lengths 2x1. It has been a subject of active study in 
combinatorics and statistical mechanics since the number of domino tilings in 
rectangular regions was found [U [2]. 

One can find an implementation of domino tilings in arrangements of full 
tatami mats on the floor of old Japanese style rooms. The floor is covered with 
1.82m X 0.91m-sized rectangular mats. Traditionally, tatami mats are arranged 
under the condition that four corners meeting at a point is avoided because the 
pattern is believed to be inauspicious. In this article, we call this additional rule 
the tatami condition. 

One can generalize the tiling by allowing I x 1 tiles or tatami half-mats. The 
tatami condition makes sense even if half-mats are present. It states no four 
corners meet regardless the type of tiles. A tiling with 2x1 and 1x1 tiles that 
satishes the tatami condition on all the vertices shall be called a tatami tiling. 

The rectangle-square tiling model without tatami condition has been stud¬ 
ied by mapping it onto the monomer-dimer models. Various interactions has 
been considered for that model, namely nearest neighbor monomer-monomer, 
monomer-dimer, dimer-dimer interactions, the one depending on the orientation 
of dimers, etc. in addition to the hard-core one. The tatami condition is realized 
as a novel nearest neighbor four-body interaction, which is worth investigation 
in its own right. As we find below, the tatami condition makes the Monte Carlo 
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sampling hard, it is a suitable model for testing techniques to improve sampling 
efficiency. 

For this model, a natural combinatorial question arises: How many tatami 
tilings are there for a given region? Erickson et. al. obtained the number for 
small rectangular regions in refs. [am [5] with the transfer matrix approach. 
Further, they studied the tilings when the number of square tiles is given. Let 
the number of tatami tilings with m squares on x £2 rectangular 
region. They obtained 




{ m • 2™ + (to + 1) ■ 

0 


(m < £,m = £ mod 2) 

(to = £) 

(otherwise). 


( 1 ) 


via a beautiful combinatorial argument [B1 [7]. The number of all the tatami 
tiling is obtained as 

Y, tiAm) = 2^-\3£ - 4) + 2. (2) 

m—0 


This expression suggests that the degrees of freedom of this statistical model 
lives on the boundary because the free energy or the logarithm of ([2]) is of order 
0{£) instead of order 0{£'^) [1]. 

The square-rectangle tiling without tatami condition has bulk degrees of free¬ 
dom: the logarithm of the number of conhgurations has order 0{£^). Therefore 
the tatami condition seems to change the nature of the system. It is desirable 
to study the interpolating systems and find the number of configurations with 
the number of tatami condition violated points given. It is interesting to ask 
whether this phenomenon occurs on other geometries. Combinatorial methods 
EllI] , however, makes full use of the strong condition: the tatami condition 
and the geometry. Therefore it is hard to apply to other situations. Methods 
applicable to general cases are desirable. 

The purpose of this work is two-fold. One is to provide an answer to an in¬ 
teresting combinatorial problem, the estimation of the number of conhgurations 
of the tatami tilings on a given region. The other is to develop a Monte Carlo 
(MC) conhguration counting technique for this specihc hardly relaxing system. 

In this article, we generate and count them stochastically by the Markov 
chain Monte Carlo (MCMC) method in contrast to combinatorial generation of 
tilings [5]. It works for arbitrary geometries as well as the simple ones for which 
exact enumeration or transfer matrix technique applies. We employ replica 
exchange Markov chain Monte Carlo [^ to count the frequency of visits to each 
energy state. Then histogram reweighting [TOl [H] is applied to obtain the 
density of states. 

This article is organized as follows. We dehne the statistical model in Sec. 2 
and describe its exact properties in Sec. 3. The method of MC simulation is 
explained Sec. 4. We show our estimation in Sec. 5 and state our conclusion in 
Sec. 6. 
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2 Model 


Let i? be a finite collection of faces (n,n + 1) x {m,m + 1) C {n,m € Z) 
of the planar square lattice. A configuration s is an assignment of domino tiles 
on R. To be precise, we put 2 x 1-sized domino tiles each of which covers two 
adjacent faces and does not overwrap with other tiles. All the faces that are not 
covered with the domino tiles shall be covered with square tiles of 1 x 1 size. 
We say that a configuration s satisfies the tatami condition on a vertex when 
less than four tiles meet there. 

We map this model onto a monomer-dimer model on the dual square lattice 
as depicted in Fig. [T] A domino tile corresponds to a dimer which is placed on 
a dual edge. A dual vertex v* is an endpoint of at most one dimer. If no dimer 
touches V* , one says that a monomer occupies v*. 
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Figure 1: (a) An example of tatami tilings with four 2x1 tiles and a 1 x 1 tile, 
(b) A tiling for which the tatami condition is broken at a vertex (indicated by 
a circle), (a)’ and (b)’ are the monomer-dimer representation of (a) and (b). 


We consider a monomer-dimer model defined with the partition function 

Z(AJt,Jd)= ^ ^ (3) 

s: monomer-dimer covering 

where /3 is the inverse temperature and the Hamiltonian is given by 


H{s- Jt, Jd) = JtNt{s) + JdlVd(s). (4) 

The quantities 

^t(s) = ^ nt{f*,s), (5) 

f* :dual face 

Nd{s) = ^ nd{v*,s) ( 6 ) 

V* :dual vertex 

are the number of dual faces where the tatami condition is broken and that of 
dimers, respectively. The coefficients Jt and Jd can be chosen arbitrarily. The 
term Nt can be regarded as a local multi-tile interaction in the original picture. 
Both can be written as sums of local quantities. 




nd{v*,s) 


if some of dual edges around /* are occupied by dimers 
otherwise 


(7) 


i if V* is occupied by a dimer 
0 otherwise 


( 8 ) 
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In this article, we always impose the hard wall boundary condition 0, that 
is, the tatami condition always holds on vertices on the boundary. 

The partition function ([3]) can be rewritten as 

Z{/3; Jt, Jd) = ^ g{N, M) ^ 

N,M 

where g{N, M) is the number of configurations with N tatami condition violated 
dual faces and M dimers. If we define 

g{N)=Y,9iN,M), (10) 

M 

the number of all the tatami tilings is given by ^(0). 

3 Exact Relations 

Though our method is applicable to arbitrary plane regions, in this article, we 
report the results of rectangular and aztec diamond regions. An aztec diamond 
of order i consists of all faces whose center {x, y) satisfy |a;| + \y\ < i. It consists 
of + 1) faces. For regions of this shape, the number of tatami tilings is 
not known exactly but , for small fixed can be counted with an algorithm 
applicable to arbitrary lattices m- 

For a fixed region, let iVmax be the largest N such that g{N) > 0. It can 
be shown that iVmax = (-^i — 1) x (^2 — 1) for x £2 rectangular region, and 
A^max = 2I'(€ — I) + 1 for the aztec diamond of order 1. For large rectangular or 
aztec diamond regions, the relation 

5(A^max) = 1 (11) 

holds. This is because iVmax is attained solely by the state s consisting only 
of monomer^. This relation will be used for normalizing the density of states 
from Monte Carlo samples in Sec. 14.21 

For some regions, we have more exact relations for g. For example, we have 

5(A'max - 1) = 2(£i - 1) + 2(£2 - 1) (12) 

for £i X £2 rectangular region. This is the number of M = I configurations for 
which the dimer is placed along a boundary. These additional relations is useful 
for the normalizing the density and for verifying the accuracy of the numerical 
calculation. 


4 Estimation Method 

We estimate the density of states g{N,M) and g{N) defined in eqs. © and 

m- 

^One could impose an alternative boundary condition, namely the open boundary condi¬ 
tion, where the compliment of R is populated with 1x1 tiles. 

^ Large regions with smooth boundaries have this property. There are a number of regions 
for which ^(A^max) is not unity but is known exactly. For example, ^(A^max = 0) = 3 for 1x3 
rectangular region. Note the hard-wall boundary condition. 
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4.1 Monte Carlo Moves 


We consider a Markov chain on the configuration space {s} which satisfies the 
detailed balance condition and has the distribution as the stationary state. 
To make a transition, we choose an edge uniformly random and apply the move 
corresponding to the local configuration as depicted in Fig. [5] Then we accept 
or reject the new configuration using Metropolis algorithm with Hamiltonian 
dH). In ref. [13], the Markov chain with the moves (a) and (b) only was used. 
We propose to add the move (c) in order to sample configurations with many 
dimers and many tatami condition violated dual faces efficiently. 



(a) 



PzSI 


(C) 



Figure 2: Set of moves, (a) If the chosen edge is the center of a dimer or two 
monomers, change the state to two monomers or a dimer, respectively, (b) If 
the edge is between a monomer and a dimer, swap them, (c) If the edge is one 
of two edges shared by two dimers, rotate the dimers 7r/2. (d) If none of (a)-(c) 
applies, do nothing. 


4.2 Estimation of g{N) 

To count the number of tatami configurations (?(0), we set Jd = 0. We also set 
Jt = +1 without loss of generality. 

We run the simulation at a series of temperatures {/3j}. Though one could 
obtain Z g{Q) by sending /3 to infinity, the histogram reweighting [TU] [2] is 
more efficient. Let be the number of visits to each energy level N = Nt 
normalized by all the visits at Pj. We note that in the limit of infinite sample 
size, /lATj converges to 

^ _ g{N) exp(-/3j7V) 

^N,j — 


zm 


due to the law of large numbers. We have 


9iN) = 


It-n, 


exp(-/3,iV)/Z(/3,) 


independent of temperatures once Z{Pj) is known. 

For finite samples, these estimates have statistical errors. In ref. a single 
estimate is proposed by composing sums in the numerator and the denominator 
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as 



( 13 ) 


where Z{fij) is determined self-consistently by 


=^9{N)exp{-^jN). 


(14) 


N 


In ref. 1 15). it is proved that this set of equations gives the maximum likelihood 
estimation of A-g{N) {N = 0,..., IVmax) up to an over-all factor A. In practice, 
we start with g{N) = 1 and Z{j3j) = 1 and solve eqs. iteratively for 

two families of unknowns g{N) and Z{l3j). The over-all factor A can be fixed 
by eq. mi)- 

It is numerically hard to do calibration at N = iVmax- High energy states 
near N = fVmax are sampled rarely due to the penalty term Jt > 0 and small 
density of states (see eqs. (ITTl) . (IT^ I. This can cause serious loss of accuracy of 
the estimate for A ■ ^(A^max)- 

One would have a large sample if Jt was negative and the tatami condition 
breaking vertices were favored. Obviously, this swaps the situations of small 
and large N regions; negative Jt makes the estimate of A ■ g{0) inaccurate. A 
cure is to measure A± ■ g{N)± for large (small) N with positive (negative) Jt 
respectively and match the estimated values as g{N„i)+ = g{N^)- at some 
intermediate value cs A^„iax/2. The problem is that there is no single way 
to choose A^ni- One could even minimize the sum of differences at several N^’s. 

Instead of matching two groups of data as above, we redefine /3 = /I • Jt = ±/3 
and consider a series of ‘temperatures’ 


(15) 


I3j = j ■ A/3, 


where j takes zero, positive and negative integer values. By running simulations 
at these ‘temperatures’, a single set of histograms is generated. This method 
is superior to the former matching method in that it is simply the maximum 
likelihood estimation and is free of ambiguity of the matching condition. The 
iterative process takes longer to converge though. 

One could say that data for ferromagnetic and anti-ferromagnetic models 
are combined to make a single family of histograms. We can apply the or¬ 
dinary reweighting equations dni), (HU) to obtain g{N) because reweighting 
works not only for thermodynamic distributions but also for general probability 
distributions [TBl [T71 ITS] . 

4.3 Replica Exchange Monte Carlo 

We speed up the simulation by the replica exchange Monte Carlo [S] as was done 
for configuration counting problem in ref. El- 

For (3 defined in eq. (HU, we consider exchanges of replica pairs at neighbor¬ 
ing temperatures, including the pairs 0j,j3j+i) = (—A/3,0), (0, -l-A/3). The pro¬ 
cedure can be regarded as a variant of Hamiltonian replica exchange method |161 
[H] where one considers replicas having different Hamiltonians. To our knowl¬ 
edge, however, there has been no application in which Hamiltonians with posi¬ 
tive and negative coefficients are exchanged. 
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4.4 Estimation of g{N, M) 

To estimate the number of configurations for the number of dimers given, we 
make both Jt and Jd vary. In this case, we set /3 = 1 and regard (Jt, Jd) as free 
parameters. We run multi-dimensional replica exchange Monte Carlo [l6l [I2] 
with replicas having the set of parameters 

Jd,k) = U ■ AJt, k ■ AJd), (16) 

where j, k takes positive and negative integer values in the spirit of (HU. 

A replica (j • A Jt, fc • AJd) exchanges its configuration with ((j ± 1) • AJt, fc • 
AJd) and (j • AJt, (fc ± 1) ■ AJd). 

Using multi-parameter reweighting [T51 118) . we obtain the multi-dimensional 

density of states g{N,M). Namely, one iteratively solve the following set of 
equations 

^ exp(-( Jt,,JV + Jd,k) ’ 

ZiJt,j,Jd,k) = ff(A^,Af)exp(-(Jt,,W +Jd,fcM)), (18) 

N,M 

which is the straightforward generalization of eas ]13ll4l The two-dimensional 
histogram hpfMjk stands for the number of visits to (W, Ad) = {N,M) state 
by the (j, fc)-th replica normalized by all its visits. 

The over-all factor is fixed using an exact relation like (HU. For example, 
5(Ai„ax, 0) = 1 can be used for rectangular regions. 

5 Results 

We perform the estimation on rectangular regions up to ii + £2 < 256 and 
aztec diamond regions up to £ < 48. To estimate the statistical error, we divide 
simulation data into 16 blocks and calculate the variation among block averages. 
Simulations have been performed on a computer with modest power. 

• CPU : Xeon E5-2640 x2 

• memory : 32GB 

• OS : Windows 7 Professional 

• compiler : Intel C-I--I- Compiler 

5.1 Estimation of g{N) and the Number of Tatami Tilings 

^( 0 ) 

As an example histogram Hnj for a square region is shown in Fig. [3] One 
confirms that configuration at all energies 0<A<(£ — 1)^ are visited thanks 
to the choice m- 

The density g{N) for a square region and that for the aztec diamond region 
are shown in Figs.|4]and[^ They are all unimodal in N for all the sizes £ studied. 
One finds the estimates for the numbers of tatami tilings as g{0) in Figs. 2] and 

El 
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Figure 3: Histogram H{E) is Hej multiplied by the sample size 5 x 10^ for 
10 X 10 square region. Each curve corresponds to j3j shown in the legend. 
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Figure 4: Estimation of the number of tilings g{E) on 10 x 10 square region. 
The over-all factor A is fixed by imposing g(81) = 1. The number of tatami 
tilings g{0) ~ 1.33 x lO"* can be read on the vertical axis. 











The estimates ( 7 ( 0 ) for square and aztec diamond regions are shown in Ta¬ 
bles [T] and 15.11 respectively. They are consistent with the exact result ([5]) for 
square regions and the result of direct enumeration for small aztec diamonds 
within statistical error. 

5.2 Size Dependence of g{N) 

In contrast to the exact result no exact formula for g{N) are available for 
> 0. We plot log g{N) against the size £ of regions in Fig. [SI It is suggested 
that \ogg{N) {N = 1, 2) are also of order 0{£) in the large i limit. 

One can understand this behavior in analogy with the low temperature ex¬ 
pansion for spin systems. A generic ground state N = 0 has 0{£) dimers along 
the boundary. If we split one of them into two monomers, we have an iV = 1 
configuration (Fig. jT]). We have more = 1 tiling that has a tatami condition 
violated point in the interior. Thus we have g{l) > O{g{0)) x 0{£). 

5.3 Dependence on Boundary Lengths 

We plot log 5 ( 0 ) against the boundary length p = 2{£i + £ 2 ) for families of 
rectangles with fixed aspect ratios and p = M for aztec diamonds in Fig. [SI We 
see that leading contribution to log (/(O) is of the linear form C-p+D. When one 
fits the data with the least square method, the estimated value of C is within 
3% of i log 2, the exact value ([2]) for square regions, for all the families. On the 
other hand, the constant term D depends on the shape and the aspect ratio and 
takes values between 2.0 and 4.5. 

5.4 Estimation of g{N, M) 

By the method in Sec. 14.41 we obtain a table of g{N, M) for each plane region 
R. As examples, data for 8 x 8 square region and the aztec diamond region of 
order 5 are drawn as heat maps in Figs. [TOl and ITT] 

The fact that an estimated value of g{N^ M) is positive implies that the 
unknown true value is in fact positive. The converse, however, is not always 
true because the simulation can fail to visit any states at {N,M). 

The density prohle g{N,M) takes a similar form for all the regions investi¬ 
gated. The density g is non-zero in a domain 

{[N,M)\N_{M) <N< N+{M)}. 

It has a maximum at a cell in the interior of the domain and is decreasing with 
the distance from that cell. Two limits N = N±{M) in {N,M) plane cross at 
{N,M) = (Afn.ax,0). 

Then the lower bound is given by A^ > A^max — 2M. This inequality can 
be shown as follows. Consider the tiling {N,M) = (Amax,0) consisting only of 
monomers. Then one introduces dimers one by one. Each dimer creates at most 
two dual faces where the tatami condition is met. This bound is saturated for 
small M for the cases investigated. 

The upper bound is given by A^ < A^max — M for small M. This comes from 
the fact that introduction of the first dimer creates, even if it is placed along 
a boundary line, at least one dual face where the tatami condition is met. For 
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Figure 5: Estimation of the number of tilings g{E) on the aztec diamond region 
at order i = 5. The over-all factor A is fixed by imposing g(41) = 1. The 
number of tatami tilings g{0) ~ 3.08 x 10^ can be read on the vertical axis. 


Table 1: Number of tatami tilings g{0) of square regions of various sizes. Esti¬ 
mation by MC and exact results are compared. Statistical error is estimated by 
dividing the whole simulation into 16 blocks. MCS represents that within one 
block. 


Size i X i Exact [3] (2^ ^(3^ — 4)-I-2) MC result MCS x # of replica 


6 

X 

6 

450 

(4.49 ±0.03) 

X 

10 ^ 

(3.0 

X 

10 ®) 

X 

48 

7 

X 

7 

1090 

(1.09 ±0.01) 

X 

10 ® 

(3.0 

X 

10 ®) 

X 

48 

8 

X 

8 

2562 

(2.57 ±0.02) 

X 

10 ® 

(4.0 

X 

10 ®) 

X 

48 

9 

X 

9 

5890 

(5.89 ±0.07) 

X 

10 ® 

(4.0 

X 

10 ®) 

X 

48 

10 

X 

10 

13314 

(1.33 ±0.02) 

X 

10 ^ 

(5.0 

X 

10 ®) 

X 

48 

11 

X 

11 

29698 

(2.96 ±0.04) 

X 

10 ^ 

(6.0 

X 

10 ®) 

X 

48 

12 

X 

12 

65538 

(6.56 ±0.06) 

X 

10 ^ 

(8.0 

X 

10 ®) 

X 

48 

15 

X 

15 

671746 

(6.70 ±0.07) 

X 

10 ® 

(2.0 

X 

10 ®) 

X 

48 

20 

X 

20 

29360130 

(2.93 ±0.04) 

X 

10 " 

(5.0 

X 

10 ®) 

X 

48 

25 

X 

25 

1191182338 

(1.19±0.01) 

X 

10 ® 

(2.0 

X 

10 ") 

X 

48 

32 

X 

32 

197568495618 

(2.00 ±0.07) 

X 

10 “ 

(1.0 

X 

10 ®) 

X 

96 

40 

X 

40 

63771674411010 

(6.38 ±0.20) 

X 

10 ®® 

(2.0 

X 

10 ®) 

X 

96 

48 

X 

48 

1.9703248369746 x lO^® 

(2.04 ±0.12) 

X 

10 ®® 

(2.0 

X 

10 ®) 

X 

96 

64 

X 

64 

1.7339939429287 X lO^i 

(1.69 ±0.16) 

X 

10 “ 

(2.0 

X 

10 ®) 

X 

96 

128 

X 

128 

6.4653649714978 x 10^° 

(1.79 ± 1.18) 

X 

10 “ 

(2.0 

X 

10 ®) 

X 

96 


10 






Table 2: Number of tatami tilings ^(0) of aztec diamond regions of order 
Estimates by MC and exact results are compared. 


Order i 

Enumeration 

MC 

MCS X # of replica 

2 

80 

(8.00 ±0.04) X 10^ 

(4.0 X 10^) X 48 

3 

392 

(3.92 ±0.03) X 10^ 

(5.0 X 10^) X 48 

4 

1200 

(1.20 ±0.01) X 10^ 

(1.0 X 10®) X 48 

5 

3080 

(3.08 ±0.03) X 10^ 

(2.0 X 10^) X 48 

6 

7312 

(7.32 ±0.06) X 10^ 

(4.0 X 10^) X 48 

7 

16712 

(1.68 ±0.02) X 10^ 

(6.0 X 10^) X 48 

8 

37424 

(3.75 ±0.05) X 10^ 

(1.0 X 10®) X 48 

10 

181392 

(1.81 ±0.02) X 10^ 

(2.0 X 10®) X 48 

12 

- 

(8.51 ±0.13) X 10^ 

(2.0 X 10®) X 48 

16 

- 

(1.78 ±0.04) X 10^ 

(2.0 X 10®) X 48 

24 

- 

(6.57 ±0.34) X 10^ 

(2.0 X 10®) X 48 

32 

- 

(2.30 ±0.17) X 1012 

(2.0 X 10®) X 48 

48 

- 

(2.17 ±0.27) X lOi’’ 

(2.0 X 10®) X 48 
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Figure 6: The number of tilings with given number of violation. The system 
size stands for the edge length t of square regions. 




Figure 7: An example of obtaining an = 1 tiling configuration from = 0 
one by splitting a dimer into two monomers. 
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Figure 8: An example of A^ = 1 configuration that cannot be obtained from 
N = 0 one by the splitting procedure in Fig. 0 



perimeter 


Figure 9: The number of tatami tilings g{0) against boundary length. The 
vertical axis stands for log 5 ( 0 ) while the horizontal axis stands for p, the quarter 
of the boundary length. Ratio 1:2, 1:3 and 2:3 mean rectangles with these aspect 
ratios. 
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aztec regions it holds only for M = 0,1 while for rectangular regions it holds 
for 0 < M < min(£i, € 2 )- 

There is a difference between the shapes of the domains for rectangles and 
those for aztec diamonds, especially near the bottom row N — 0 which cor¬ 
responds to tatami tilings. For square region case, the result for the row is 
consistent with ©• For aztec diamond region case, the data for the row is 
consistent with the result of direct enumeration; at order 5, the tatami tilings 
with less than two monomers are absent. 


6 Conclusion 

We have studied a monomer-dimer model with an interaction term which en¬ 
forces the tatami condition. We have been able to estimate the number of 
the tatami tilings and find the size dependence. Moreover, by two-dimensional 
replica exchange Monte Carlo and multi-parameter histogram reweighting, we 
have been able to estimate the number of tilings with specified numbers of 
dimers and tatami condition violated dual faces. To this end, we have proposed 
a Monte Carlo method to calculate the density of states of combinatorial models 
by combining ferromagnetic and anti-ferromagnetic models. Efficiency of this 
method when applied to other combinatorial problem is left for future study. 
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Figure 10: Heat map of estimated logiQ g{N,M) for the 8x8 square region. 
The vertical and the horizontal axes represent N and M, respectively. For white 
cells, the numbers of tilings are estimated to be zero. 
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Figure 11: Heat map of estimated log^g g{N, M) for the aztec diamond region of 
order 5. The vertical and the horizontal axes represent N and M, respectively. 
For white cells, the numbers of tilings are estimated to be zero. 
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